---
title: "Lab 1 Code"
author: "Fabio Setti"
date: "Spring-2026"
output: html_document
---

## Today's Packages 🤗

```{r}
install.packages("rio")
install.packages("psych")
install.packages("tidyverse")
install.packages("car")
```


```{r}
library(car)
library(psych)
library(tidyverse)
```


```{r}
dat <- rio::import("https://fabio-setti.netlify.app/data/NY_Temp.txt")
```

## Describing the Data Numerically


```{r}
str(dat)
```


```{r}
# the `trim = .05` argument calculates means without the bottom/top 2.5% of the data
# add IQR = TRUE to get the interquartile range
psych::describe(dat[,-1], trim = .05)
```


## Describing the Data Graphically


```{r}
ggplot() 
```


## Adding Coordinates


```{r}
ggplot(dat,  
 aes(x = Case, y = Temp)) 
```


## A Basic Plot


```{r}
ggplot(dat,  
 aes(x = Case, y = Temp)) +
 geom_point()
```

## Change the Theme

```{r}
ggplot(dat,  
 aes(x = Case, y = Temp)) +
 geom_point() +
 theme_classic()
```


```{r}
# set theme globally
theme_set(theme_classic())
```

## Add Linear Regression Line

```{r}

ggplot(dat,  
 aes(x = Case, y = Temp)) +
 geom_point() +
 geom_smooth(method = "lm")
```


## Add a Loess Line


```{r}
ggplot(dat,  
 aes(x = Case, y = Temp)) +
 geom_point() +
 geom_smooth(method = "lm", se = FALSE) +
 geom_smooth(method = "loess", se = FALSE,
             col = "red")
```


## 1D Scatterplot 🤔 


```{r}
ggplot(dat,  
       aes(x = Temp, y = 0)) +
       # shape = n can be used to change the shape of the dots. there are 25 shapes (1 to 25)
       geom_point(shape = 1) +
       theme_classic()
```


## Histograms


```{r}
ggplot(dat,  
       # note that we only need to give X
       aes(x = Temp)) +
       geom_histogram()
```


## Histogram: Adjusting Bins


```{r}
ggplot(dat,  
 aes(x = Temp)) +
 geom_histogram(bins = 15,
                color = "black",
                linewidth = .8,
                fill = "#3492eb") 
```


## Box-Plots


```{r}
ggplot(dat,
       aes(y = Temp)) +
       geom_boxplot()
```


## Kernel Density plots 

```{r}
ggplot(dat,
       aes(x = Temp)) +
       geom_density() +
       xlim(45, 110)
```


## Change Kernel and Smoothing Parameter


```{r}
ggplot(dat,
 aes(x = Temp)) +
 geom_density(kernel = "triangular",
              adjust = 0.5) +
 xlim(45, 110)
```


## QQplots

```{r}
ggplot(dat,
       aes(sample = Temp)) + 
       stat_qq() + 
       stat_qq_line()
```


# One-Predictor Regression

## Our Variables

```{r}
ggplot(dat,  
 aes(x = Wind, y = Temp)) +
 geom_point() +
 geom_smooth(method = "lm", 
             se = FALSE)
```


```{r}
cor(dat$Wind, dat$Temp)
```

## Run a Regression with `lm()`


```{r}
reg <- lm(Temp ~ Wind, data = dat)
```


```{r}
class(reg)
```


## What's inside the `lm()` output? 🧐


```{r}
names(reg)
```


## The `summary()` Function 


```{r}
sum_reg <- summary(reg) # normally I would not save the summary as an object, but I want to show the contents later
sum_reg
```


```{r}
class(sum_reg)
```


## What's inside `summary(lm)`? 🧐


```{r}
names(sum_reg)
```

## Making Predictions The "right way"


```{r}
reg$coef
```


```{r}

reg$coef[1] + reg$coef[2]*13 #  I prefer `as.numeric(coef(reg)[1] + coef(reg)[2]*13)` which makes sure to unnamed the output; this happens because reg$coef if a named vector, so the names remain unless we use `as.numeric()`
```


```{r}
# first we create new data that we want to get predictions on
pred_dat <- data.frame("Wind" = seq(10, 20, by = 0.2))
```


```{r}
pred_dat$y_pred <- predict(reg, pred_dat)
str(pred_dat)
```

## Visualizing Predictions


```{r}
ggplot(dat,
 aes(x = Wind, y = Temp)) +
 geom_point(shape = 1) +
 geom_point(aes(y = reg$fitted), 
                col = "red")
```

## Visualizing Residuals

```{r}
ggplot(dat,
 aes(x = Wind, y = Temp)) +
 geom_point(shape = 1) +
 geom_point(aes(y = reg$fitted), 
                col = "red") +
 geom_segment(y = reg$fitted,
              x = dat$Wind,
              xend = dat$Wind,
              yend = dat$Temp,
              linetype = 2,
              alpha = 0.5)
```

## Normality of Residuals

```{r}
car::qqPlot(reg)
```

## Plotting Residuals Against Predictors


```{r}
ggplot(dat, 
  aes(x = Wind, y = reg$resid)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_smooth(method = "loess",
              se = FALSE)
```

# Appendix

## Understanding QQplots


```{r}
# quantile values of .1, .5., .8 for normal distribution with mean = 0 and sd = 1
# this will results in Z-scores that correspond to the 
# 10th, 50th, and 80th percentile

qnorm(c(.1, .5, .8), mean = 0, sd = 1)
```


```{r}
# The quantile function calculates the observed percentiles fo a given variable.
# by supplying c(.1, .5, .8) we want the 10th, 50th, and 80th percentile of Temp

quantile(dat$Temp, c(.1, .5, .8))
```

## Creating QQplot Data


```{r}
# first, it is useful to define all the quantiles that we want as a vector
quantiles <- seq(from = .01, to = .99, by = .01)  
str(quantiles)

# calculate quantiles for normal distribution and real data
Xaxis <- qnorm(quantiles)
Yaxis <- quantile(dat$Temp, quantiles)

# Note that you can use ANY mean and SD value for the qnorm() function.
# Here, I use the mean and SD of the Temp variable
Xaxis_2 <- qnorm(quantiles, mean = mean(dat$Temp), sd = sd(dat$Temp))

# create a data.frame (ggplot like data.frame objects)
QQdata <- data.frame("Xaxis" = Xaxis,
                     "Xaxis_2" = Xaxis_2,
                     "Yaxis" = Yaxis)
head(QQdata)
```

## QQplot The Data

```{r}
ggplot(QQdata, 
       aes(x = Xaxis, y = Yaxis)) +
      geom_point()
```


## A more intuitive "QQplot"?


```{r}
ggplot(QQdata, 
       aes(x = Xaxis_2, y = Yaxis)) +
       geom_point() +
       # equivalent to stat_qq_line()
       geom_abline(intercept = 0, slope = 1)
```